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ABSTRACT 

A significant fraction of white dwarfs show metal lines indicative of pollution with planetary 
material but the accretion process remains poorly understood. The main aim of this paper 
is to produce a road-map illustrating several potential routes for white dwarf pollution and 
to link these paths to observational outcomes. Our proposed main road begins with the tidal 
disruption of a scattered asteroid and the formation of a highly eccentric tidal disc with a 
wide range of fragment sizes. Accretion of these fragments by Poynting-Robertson (PR) drag 
alone is too slow to explain the observed rates. Instead, in the second stage, several processes 
including differential apsidal precession cause high-velocity collisions between the eccentric 
fragments. Large asteroids produce more fragments when they disrupt, causing rapid grind- 
down and generating short and intense bursts of dust production, whereas smaller asteroids 
grind down over longer periods of time. In the final stage, the collisionally produced dust 
circularises and accretes onto the white dwarf via drag forces. We show that optically thin dust 
accretion by PR drag produces large infrared (IR) excesses when the accretion rate exceeds 
107 g/s. We hypothesise that around white dwarfs accreting at a high rate, but with no detected 
infrared excess, dust circularisation requires enhanced drag - for instance due to the presence 
of gas near the disc's pericentre. 
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1 INTRODUCTION 


Over the last decade, it has become clear that planetary systems 
around solar mass stars are the norm rather than the exception, with 
estimated occurrence rates around unity (Berta-Thompson et al. 
2015; Dressing & Charbonneau 2015; Mulders et al. 2019; Zhu & 
Dong 2021). When these stars exit the main-sequence, subsequent 
stellar mass-loss leads to an expansion of the surrounding orbits, 
protecting material outside a few AU from the increased stellar flux 
(Veras et al. 2011; Veras & Tout 2012; Veras 2016), while any ob- 
jects within are engulfed (Mustill & Villaver 2012; Villaver et al. 
2014). The enduring presence of planetary material around white 
dwarfs is inferred from the fraction of 2596 - 5696 of systems that 
show metal lines in their spectra (Zuckerman et al. 2003, 2010; 
Koester et al. 2014; Wilson et al. 2019). So far, the presence of 21 
different heavy elements has been detected in white dwarf photo- 
spheres, with 19 of these found in the single system GD 362 (Zuck- 
erman et al. 2007; Xu et al. 2013, 2017; Melis & Dufour 2017). 
Detailed analysis of these polluted white dwarfs has opened a 
new channel for the investigation of exoplanetary systems. In par- 
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ticular, the detection of multiple heavy elements provides an op- 
portunity to directly study the composition of planetary material 
outside the Solar System. In most cases, both the abundances (Xu 
et al. 2014, 2019; Wilson et al. 2015) and oxygen fugacity (Doyle 
et al. 2020, 2021) are found to match bulk Earth to zeroth order, 
although core-rich, mantle-rich and even crust-dominated objects 
have also been observed (Hollands et al. 2017; Swan et al. 2019b; 
Hollands et al. 2021). These systems have been interpreted as evi- 
dence for early planetsimal formation followed by collisional pro- 
cessing (Harrison et al. 2018, 2021) and indicate that the major- 
ity of exo-planetesimals are differentiated (Bonsor et al. 2020). In 
future models, accounting for pressure-sensitivity of different ele- 
ments will allow the origin of the accreted objects to be determined 
as well (Buchan et al. in prep). 


However, while the study of polluted white dwarfs has blos- 
somed into a field of its own, the processes via which planetary 
material is accreted remain poorly understood. We know that the 
initial trigger for pollution is likely to be stellar mass loss, which 
widens planetary orbits and strengthens the interactions between 
planets. This can destabilise tightly-packed planetary systems, even 
if they were previously stable (Debes & Sigurdsson 2002; Maldon- 
ado et al. 2020, 2021), while systems with more space between the 
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Figure 1. Schematic diagram of the main road to white dwarf pollution examined in this work. In the precursor to its pollution, the star sheds its outer layers 
during post-main sequence evolution (dashed border). This widens and destabilises the orbits of surrounding bodies and causes some asteroids to be scattered 
towards the star through interactions with nearby planets. In the first stage of accretion, asteroids that cross the Roche radius are tidally disrupted and form 
eccentric structures with an orbital spread, which we refer to as (eccentric) tidal discs. The orbits of the surviving fragments are perturbed via various processes, 
including differential precession, causing high-velocity collisions and grind-down within the eccentric tidal disc (stage II). The resulting dust then circularises 


and accretes onto the star via various drag forces (stage III). 


planets likely survive intact (Duncan & Lissauer 1998; Veras et al. 
2016). If the planetary system contains an asteroid belt, a single 
massive planet orbiting interior to the belt can scatter large num- 
bers of asteroids within its expanding chaotic zone (Bonsor et al. 
2011; Mustill et al. 2018). An outer planet can be similarly effec- 
tive and scatter asteroids around expanding interior mean motion 
resonances (Debes et al. 2012; Antoniadou & Veras 2016). Aster- 
oids or planets that pass within the Roche radius tidally disrupt into 
highly eccentric discs whose shapes range from narrow if the as- 
teroid was small (Veras et al. 2014, 2021; Nixon et al. 2020), to 
wider and even partially unbound if the object was terrestrial-sized 
(Malamud & Perets 2020a,b). 


From this point on, the evolution of the fragments remains 
more obscured. As an end-point of their evolution, observations 
show that a minority of polluted white dwarfs are surrounded by 
dust (Rocchetto et al. 2015; Farihi 2016; Wilson et al. 2019), whose 
emission typically varies within several years at mid-infrared (Far- 
ihi et al. 2018; Swan et al. 2019a, 2020), but less commonly at 
near-infrared (Rogers et al. 2020). Some of these systems also show 
evidence of on-going gas production (Gánsicke et al. 2006, 2007, 
2008; Manser et al. 2020) and they occasionally contain larger, 
transiting bodies (Vanderburg et al. 2015; Manser et al. 2019; Van- 
derbosch et al. 2020, 2021; Guidry et al. 2021). While no complete 
description of the accretion process currently exists, the general 
hypothesis is that small fragments are circularised by Poynting- 
Robertson (PR) drag (Rafikov 2011b; Veras et al. 2015a,b), while 
larger fragments require a prior phase of collisional grind-down 
(Jura 2003; Jura et al. 2007; Wyatt et al. 2011; Li et al. 2021). Other 
mechanisms that induce orbital changes after disruption are the 
Yarkovski force (Veras et al. 2015a,b; Malamud & Perets 2020b), 
potential interactions with pre-existing material around the star (Gr- 
ishin & Veras 2019; Malamud et al. 2021) and magnetic Alfvén- 
wave drag (Zhang et al. 2021). 

In this paper, we systematically investigate how planetary ma- 
terial accretes onto white dwarfs. We eventually produce a road- 
map illustrating several potential routes for white dwarf pollution 
with links to observational outcomes (Fig. 18). The main path to 
accretion that we examined (see Fig. 1) begins with the tidal dis- 
ruption of a scattered asteroid to form a highly eccentric tidal disc 
(Sect. 2). We evaluate its morphology and constrain the upper and 
lower bounds of their fragment sizes due to radiative and tidal 


forces. Then, we present a short intermezzo where we consider the 
merits and limitations of a simple collision-less evolution model via 
PR drag (Sect. 3), which we find cannot drive sufficiently rapid ac- 
cretion, even under the most optimistic assumptions. We continue 
our main road to pollution in Sect. 4, where we discuss how var- 
ious processes induce high-velocity collisions between fragments. 
These collisions take place while the fragments still travel along the 
highly eccentric orbits (e > 0.999) on which they are released. This 
notion is fundamentally different from previous models that calcu- 
lated collisions between fragments that were already supposed to 
have circularised (Kenyon & Bromley 2017a,b; Swan et al. 2021). 
In Sect. 5, we present a simple but quantitative calculation of ec- 
centric collisional grind-down based on the fragment’s differential 
rates of apsidal precession, a process that likely induces collisions 
in all tidal discs on this scale. We then enter the third stage of our 
accretion scenario where we model the geometry and emission of 
the dust that is produced (Sect. 6). We discuss variations of this 
model as well as alternative paths to white dwarf pollution in Sect. 
7, including observational outcomes when they are sufficiently well 
understood. Finally, we conclude our work in Sect. 8 


2 STAGE I: FROM ASTEROID TO TIDAL DEBRIS DISC 
2.1 The tidal disruption criterion 


The outer edge of the disruption zone is set by the distance where 
an asteroid's internal strength and self-gravity are overcome by 
tidal forces. The details of this process depend on the asteroid's 
shape and composition, as well as on its path and potential rota- 
tion (Dobrovolskis 1982, 1990; Davidsson 1999; Davidsson 2001). 
We content ourselves here by by considering an idealised case of 
breakup by tensile failure, likely the most common type of tidal 
disruption for solid bodies. We adopt a similar approach as Bear & 
Soker (2015) and Brown et al. (2017) and identify the breakup cri- 
terion for a spherical, non-rotating asteroid of size Rast and density 
Past as the point where the summed forces from the tensile strength 
(Fs) and self-gravity (zc) first fail to compensate the tidal force 
induced by the gravitational gradient (Fr): 


Fs + Fog + Fr c —SR2. 


GM? n SE E 0, (1) 
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where G is the gravitational constant, is r is the distance to the 
white dwarf and S is the material’s tensile strength. The masses of 
the white dwarf and the asteroid are indicated by Mwp and Mast, 
respectively. In the gravity-dominated regime (|Fsg| > |Fs|), the 
distance at breakup is mostly independent of the asteroid’s size and 
occurs at the classical Roche radius (Davidsson 1999; Bear & Soker 
2013; Veras et al. 2014; Malamud & Perets 2020a,b, e.g.): 


1 

2 3 

reste = ( we Rwp. Q) 
ast 


To quantify Eq. 2, it is necessary to specify the white dwarf density 


Pwo. Neglecting the slight temperature dependence, this relation- 
ship can be approximated by (Nauenberg 1972): 
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The Roche radius amounts to roughly 1 Re for a 0.6 Mo white 
dwarf. If they are monolithic (as opposed rubble-pile aggregates, 
see below), smaller asteroids can survive a certain distance within 
the Roche radius until the extra barrier of their internal strength 
is overcome. Accounting for this, Eq. 1 can be solved to yield a 
maximum object size (Rmax) that can survive at a distance r from a 
WD: 


Nie 


Pet toes sd] 


(4) 


max 


Chondrite asteroid samples indicate that the upper end of ten- 
sile strengths is around 0.1-10 MPa (Scheeres et al. 2015; Ve- 
ras & Scheeres 2020; Pohl & Britt 2020) but it is understood to 
vary by orders of magnitude depending on an asteroid's formation 
history and composition. For instance, estimates from modeling 
of cometary material composed of ice-coated interstellar silicate 
grains indicate sub-kP strengths (Greenberg et al. 1995; Davidsson 
1999; Gundlach & Blum 2016). A formation via the gentle sticking 
of constituent particles can even result in so called rubble piles with 
near-zero effective strength. Such an object is for instance invoked 
to explain the rapid breakup of Shoemaker-Levy 9 in Jupiter's outer 
envelope (Asphaug & Benz 1994). To visualise these differences, 
we plot the maximum intact object size as a function of distance 
for a range of tensile strengths in Fig. 2. The figure illustrates the 
dichotomy that arises based on the asteroid's size. Large (> 100 
km) asteroids always fragment at the Roche radius regardless of 
their tensile strength. In contrast, because the tidal force scales as 
Fry e R^ and Fs œ RÈ, smaller km-sized granite rocks can reach 
as close as a few percent of the Roche radius before they are torn 
apart. 


2.2 Tidal debris disc morphology 


As the main body begins to break up, its fragments stray from the 
initial orbit and spread out over a range of energies. Malamud & 
Perets (2020a,b) simulated this process for terrestrial-sized bodies 
and showed how the breakup typically proceeds in stages. Pericen- 
tre passages that are close to the Roche radius typically result in 
partial disruption and require several orbits in order to completely 
destroy the object. Each passage can add spin to the surviving ob- 
ject and progressively weaken it. The fragments that break off begin 
to form a stream that can gravitationally collapse perpendicular to 
its direction of motion. After several orbits, this produces a fully 
formed tidal disc of interlaced elliptical annuli. 
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Figure 2. The maximum sizes of objects that are safe from tidal disruption 
at a distance r from a 0.6 Mo white dwarf, plotted for three different ten- 
sile strengths. Large asteroids (2100 km) always disrupt close to the Roche 
radius of the star, whereas smaller fragments («€ km) can survive deep peri- 
centre passages (r < rgocne) due to their material strength. 


While the details of the breakup process lead to some varia- 
tion in the tidal disc, its basic morphological features agree with the 
simple impulse approximation method. In this framework, the dis- 
ruption is assumed to occur instantaneously at a location rg from 
the star, which for simplicity we take at the pericentre of the ob- 
ject's orbit. As the object disintegrates, its fragments are no longer 
guided by the centre of mass and continue on orbits corresponding 
to their energy at the point of separation. The ones facing the white 
dwarf are more gravitationally bound than their velocity warrants 
and move to a tighter orbit, while those on the other side of the 
asteroid migrate away from the white dwarf. Neglecting the small 
effect of binding energy in asteroid-sized objects, their specific en- 
ergy (£j) can be expressed as a sum of the kinetic (e į) and potential 
(ug) parts: 


& = &ituGi (Sa) 


1 2 1 GM: 
= =GMwp ( ) GN (5b) 
2 mp ao ri 


where rg is the breakup distance of the asteroid's centre to the white 
dwarf and ri = rg +xRagt is the distance corresponding the individ- 
ual fragments (with —1 < x < 1) and aj is the asteroid's semi-major 
axis. Correspondingly, the fragment's new semi-major axes (aj) are 
spread along 


NU AM (6a) 
. —1 
-a (122,875) (6b) 
TB 


with eccentricities (e;) equal to 


“= (1-2), (7) 
di 


The asteroid moves parallel to the plane of motion of its centre-of- 
mass prior to its breakup. This means that fragments get imparted 
an inclination (ij) depending on their vertical position that varies 
between 0 < i < Rast/r;. With this, the tidal debris discs have an 
approximate height of 2Rast at pericentre, which grows to many 
times the size of the body at apocentre, where it can be estimated at 
Hapo ~ 4Rastao /rp. In our subsequent modeling, we take the incli- 
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Figure 3. Orbital width (6a) that contains 90 % of the bound fragments 
around a 0.6 Mo white dwarf. Small asteroids on tight orbits experience 
non-dispersive breakup and form a non-dispersive stream. The width of the 
tidal disc increases with the progenitor size and its orbital separation. Very 
large objects like planets with R > Ra experience bimodal breakup with 
half of their mass being ejected from the system. 


nation constant as a function of time, although the vertical evolution 
of these tidal discs currently remains poorly understood. 
Depending on the size and semi-major axis of the asteroid 
prior to breakup, the tidal discs described as by Eqs. 6a - 7 form 
with various shapes. Asteroids that originate from wide orbits are 
only loosely bound to the star and when their size exceeds the 


threshold of Rerit = -n (Malamud & Perets 2020a), some of 
its fragments become unbound after breakup (a; < 0), with up to 
half of the material being expelled from the system in the most 
extreme case. These highly dispersive tidal disruption events have 
previously been linked to planetesimal seeding of other systems 
(Rafikov 2018). The fragments that remain on bound orbits spread 
out over a range in semi-major axes, approximately distributed 
evenly in the energy range of Eq. 5b (Malamud & Perets 2020a). 
Since the orbital energy scales as a^ !, most bound fragments be- 
come clustered on tight orbits. Therefore, we can use the inner (ain) 
and outer (aout) bound fragments to define an effective width da of 
the tidal disc as 


ac X4in (Aout — din) 
Xin + (1 — X)dou ` 


(8) 


where y is the mass fraction of bound fragments included in the 
width of the disc. We plot this effective disc width for a range of 
semi-major axes and asteroid sizes in in Fig. 3, which shows the 
clear dichotomy between the non-dispersive and bimodal regimes. 
We also visualise three examples of the different regimes in Fig. 
4. 'The top panel (a) indicates the tidal disc that forms when a 
small asteroid (1 km from 3 AU) disrupts, similar to the models 
by Debes et al. (2012), Veras et al. (2014), and Nixon et al. (2020). 
As described by these authors, the result is a completely bound 
but spaghettified orbital structure. If the size of the asteroid and its 
semi-major axis are increased (panel b), the orbital band broadens 
until ĝa >> ag and almost half of the original material is ejected 
(panel c). In this bimodal regime, any further increase of the im- 
pactor size or semi-major axis begins to increase the concentration 
of orbits closer to the star and effectively reduces the width of the 
tidal disc (see top part of Fig. 3). 
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Figure 4. Post-breakup orbits after a strengthless asteroid is tidally dis- 
rupted near a 0.6 Mo white dwarf. Panel (a) indicates the tidal disc that 
forms when a small (1 km from 3 AU) asteroid disrupts, yielding a thin or- 
bital spread. Panel (b) corresponds to a 50 km asteroid from 3 AU, which 
produces a wider - but still completely bound disc. When both the aster- 
oid size and semi-major axis are increased to 500 km and 30 AU (c), the 
outcome is a bimodal disruption event, with nearly half of the material be- 
coming unbound from the white dwarf. 


2.3 Debris size distribution 


Immediately after the main body breaks up, its fragments become 
exposed to the same tidal force that broke up their progenitor. These 
fragments are smaller than the original body and thus require a 
lower tensile strength to resist further breakup and initially remain 
stable. However, if the asteroid was scattered onto an orbit with its 
pericentre closer to the star than its initial original breakup distance, 
the fragments experience an increased peak tidal force and may dis- 
rupt again. In this simple picture, the maximum fragment size Rmax 
follows from Eq. 4 with r = reen and can be identified from Fig. 2. 
Because asteroids are most likely to be scattered to orbits with peri- 
centres near the edge of the Roche radius (Veras et al. 2021), it may 
be expected that the largest fragment is typically similar to the size 
of the asteroid itself. It is not necessary that this happens in prac- 
tice, however, as tighter orbits or rubble pile asteroids composed 
of smaller constituent particles may not lead to any large surviv- 
ing fragments. In addition, the fragments can be spun up, making 
them easier to break up (Malamud & Perets 2020b). In our model, 
we take Rmax as a free parameter due to the large associated un- 
certainty, with a baseline value of 1 km, corresponding to the size 
where fragments with a 1 kPa strength survive near the Roche ra- 
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dius. We assume the same density for the fragments as we do for 
the asteroid. 

On the other end of the distribution, the lower limit for bound 
fragments corresponds either to scale of the smallest dust grains 
that splinter off during breakup or to the blow-out size where radi- 
ation pressure from the stellar luminosity Lwp pushes the dust out 
of the system. This blow-out criterion provides an absolute lower 
limit to the fragment sizes and can be estimated from the ratio of 
radially oriented forces (Burns et al. 1979): 


Fad = 0.0013 R zi Pfrag CS Lwp 
FG um 2.7 g/cm? 0.01 Lo 


Mwp m «Qo 
CA ( 1 ) S 


where < Q > is the radiation pressure coefficient, averaged over 
the stellar emission. While white dwarfs are not luminous enough 
to blow out micron-sized grains on circular orbits, the near-unity 
eccentricities of the fragments after the tidal disruption make these 
grains susceptible. We can estimate a typical blow-out size by tak- 
ing the orbital parameters of the original asteroid with a breakup 
point at its pericentre. The smallest fragments are placed on un- 
bound orbits when when Fj4a/ Fa > 0.5(1— e), which, when com- 
bined with Eq. 7, reduces to: 


-1 
= a0 TB Cep 
Rblow = 1.51 Um (Xu) EA (a i.) 


Mwp m Ptrag zs (10) 
0.6 Mo 2.7 g/cm? 


assuming geometric scattering (< Q >= 1), which is valid for 
grains larger than the reduced peak stellar wavelength Jet /2% ~ 
0.02 — 0.1 um, depending on the stellar temperature. We plot the 
blow-out size across a range of white dwarf temperatures for three 
asteroid semi-major axes in Fig. 5. The figure indicates that the 
sizes typically range between 1-100 um depending on the white 
dwarf temperature and the asteroid's orbit. We note, however, that if 
they are produced during the tidal disruption, a population of grains 
much smaller than Rpjow can still remain remain in the system due 
to their reduced interaction with light at stellar wavelengths. We 
do not consider these tiny grains here because their properties also 
make them resistant to PR drag. 

The distribution of fragments between Rmin and Rmax is deter- 
mined by the fracture lines, the amount of sequential breakups and 
mergers, as well as the short collisional phase that follows (Mala- 
mud & Perets 2020a) and currently remains poorly constrained. We 
therefore opt to insert a truncated power-law for the fragment sizes: 


— CR, (11) 


where œ is the scaling factor of the size distribution and the con- 
stant C is set by mass conservation. Because Rmax >> Rmin and as- 
suming that o < 4, it can also be written as 


dM (4— 0) foound Mast ( R ) 
R ? 


(12) 


dR Rmax 


max 

where fpound is the mass fraction of post-breakup fragments that 
is bound to the white dwarf, a factor that follows from Eq. 6b. In 
the case of a scale-independent collisional cascade, the power law 
is characterised by œ = 3.5 and the mass is dominated by large 
fragments while the smaller particles dominate the cross section 
(Dohnanyi 1969; Tanaka et al. 1996; Wyatt et al. 2007, 2011). We 
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Figure 5. The smallest dust grains (Ryjoy) that interact with the stellar light 
(< Q >= 1) and can remain bound despite the star's radiation pressure. 
The figure assumes the disruption of a strengthless asteroid at the Roche 
radius of a 0.6 Mo white dwarf (from Eq. 10). The three lines correspond 
to different asteroid semi-major axes and the tidal disruption is assumed to 
occur at the orbit's pericentre. 


take this as our default value in our subsequently presented calcula- 
tions but note that simulations of collisional cascades that account 
for scale-dependent effects suggest a slightly lower value of alpha. 
If the mass is instead more evenly distributed over the size bins, the 
value of & is closer to 3. 


3 INTERMEZZO: COLLISION-LESS EVOLUTION VIA 
PR DRAG 


Before we consider collisions between larger fragments, we first 
evaluate the potential of perhaps the simplest scenario, where frag- 
ments accrete via PR drag alone. We derive accretion rates as a 
function of the bounded size distribution and point out the limita- 
tions of this simple scenario. 


3.1 PR contraction timescale and accretion rate 


The main contribution from PR drag on a highly eccentric orbit 
Occurs near the pericentre where the stellar flux is highest. In this 
calculation, we will make the a-priori assumption that the tidal disc 
is optically thin, which we later evaluate in Sect. 6. The averaged 
orbital equations of motion are described by Veras et al. (2015a,b): 


da. 3 « Q» LwpQ 32) 


« ; (13a) 
dt l6zpgasRac?(1— e?) 3 
d IS<Q>L 
Zoos Heng — (13b) 
dt 32z pg, Ra?c?(1— e2)2 
where we again substitute < Q >= 1 and use the simple relation 
" —1.18 
Lwp œ~ 3.26 Lo (va $ ww) (14) 
Myr 


to relate the white dwarf’s luminosity to its age fwp based on Mes- 
tel theory (Mestel 1952), with the same parameters as in Bonsor & 
Wyatt (2010). This prescription is valid up to 9 Gyr when the white 
dwarf undergoes crystallization and the cooling slows dramatically 
(Althaus et al. 2010). While more detailed cooling codes are avail- 
able (Salaris et al. 2013), Mestel’s relation captures the essential 
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Figure 6. Collision-less mass accretion rates of asteroid fragments via PR 
drag onto a 0.6 Mo white dwarf. The top panel (a) shows how the accretion 
rate varies with the size distribution (different colors) and orbital separation 
(line style). The lower panel (b) indicates what fragments reach the star via 
PR drag after a certain time (facc). Accretion starts with the smallest bound 
fragments Rplow, as bigger fragments take longer to circularise. The curves 
flatten over time due to the declining luminosity of the cooling star (Eq. 
14), meaning that fragments larger than 1 — 10 cm cannot be accreted via 
PR drag alone. 


cooling trend for the first few Gyrs, which is sufficient for our pur- 
poses here. 

The equations of motion (Eqs. 13a-14) are coupled and gen- 
erally need to be solved numerically to obtain the accretion time as 
a function of fragment size tac, (R). Because a white dwarf's lumi- 
nosity remains approximately constant for a similar period of time 
as its age, a fragment's accretion time in the window tace < twp is 
proportional to its size. For a size distribution 3 < a < 4, this leads 
to an accretion rate of: 


* — dM dtace z 
Mpr = ll dR ) = 
— (4-2) fima Masi. (15b) 
facc,max 


where facc,max is the accretion time of the largest fragment. We plot 
the PR accretion rates (accounting for stellar cooling) for 100 km 
asteroids that originate from 3 and 10 AU in Fig. 6, assuming three 
different values of a. In agreement with Eq. 15b, the accretion rate 
declines as a function of time unless @ < 3. The smallest fragments, 
therefore, typically determine the peak accretion rate in a collision- 
less scenario, even if most of the mass is contained in large frag- 
ments. The steeper the fragment size distribution, the higher the 
peak accretion rate and the steeper its decline as a function of time. 
Furthermore, we find that the accretion rate is only marginally de- 
pendent on the orbital parameters of the fragments. For a given 
value of a, the PR accretion rate can also be be written as a simple 
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Figure 7. Peak collision-less accretion rates via PR drag onto a warm (104 
K) 0.6 Mo white dwarf as a function of asteroid radius and fragment size 
distribution. The solid and (dashed) lines correspond to asteroid semi-major 
axes of 3 and (10 AU). Explaining observed implied accretion rates beyond 
10? g/s (Farihi et al. 2012) by PR drag alone requires large asteroids (Rast > 
500 km) to break up into a steep fragment size distribution (o > 3.6). 


scaling function. For the standard case of œ = 3.5, this is: 
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If the fragment size distribution is instead described by the steeper 
value of a = 3.8, the PR rate becomes: 
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(16) 
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3.2 Peak accretion rates by collision-less PR drag 


The peak accretion rate by PR drag occurs soon after the asteroid 
disrupts and the smallest fragments with size Rplow (Eq. 10) begin 
to reach the white dwarf. We evaluate these peak rates in Fig. 7 
for a range of asteroid sizes and size distributions. We take a warm 
white dwarf with temperature Twp = 104 K, which is characterised 
by both rapid PR-circularisation compared to cooler stars but also 
a larger blow-out size (Eq. 10). Of these two, the luminosity effect 
is more important so the plotted rates can be seen as upper values 
that decrease slightly for cooler stars. Fig. 7 shows that collision- 
less accretion via PR drag can at least briefly supply the lowest 
detectable accretion rates around 10° g/s in the standard case of 
a = 3.5. The higher values of observed inferred accretion rates be- 
yond 10? g/s (Farihi et al. 2012) are more difficult to explain with 
PR drag alone and require large asteroids (Rast > 500 km) to break 
into sufficiently small particles (œ > 3.6). In addition, the problem 
of explaining the high accretion rates by PR drag alone is exacer- 
bated for cooler stars with lower luminosities. 

In any case, PR drag alone is unable to accrete fragments 
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above ~ 10 cm before the star cools down (see Fig. 6). Most of 
the fragment mass is contained in much larger bodies unless the 
size distribution is incredibly steep (œ > 4). If we consider the 
default value of œ = 3.5 and a maximum fragment size of 1 km, 
the total mass fraction that can be accreted via PR drag is only 
(10 cm/Rmax) 1/2 _ 0.01. Because of this accretion inefficiency, the 
model is in conflict with the higher observed average accretion rates 
for some DBZ white dwarfs with longer sinking timescales (Farihi 
et al. 2012) and an additional process is required to explain the ob- 
served accretion rates. We will examine how collisions can play a 
role in increasing accretion efficiency in the next sections. 


4 STAGE IIa: COLLISIONS INDUCED BY ORBITAL 
PERTURBATIONS 


We propose the collisional grind-down of larger fragments into dust 
as a solution to the difficultly in obtaining high accretion rates. In 
this section, we discuss several processes that naturally lead to the 
required collisions between fragments soon after the eccentric tidal 
discs are formed. 


4.1 Differential geodetic precession 


Collisions between fragments only occur when their relative orbits 
are altered sufficiently from their initial trajectories that they begin 
to overlap. One way in which these changes could be induced, is 
from the fact that orbits do not precisely follow Newtonian tracks. 
Instead, their pericentres precess over many orbits according to GR 
(to lowest order) at a rate (GR) of (i.e. Ragozzine & Wolf 2009; 
Venkatraman Krishnan et al. 2020): 


1 
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which can lead to the build-up of apsidal differences between frag- 
ment orbits over time. In the context of white dwarf pollution, dif- 
ferential apsidal precession was first mentioned by Debes et al. 
(2012) and further examined by Veras et al. (2014), who suggested 
it can be triggered by orbital differences induced by PR drag. Be- 
cause the precession rate is highest for the most eccentric orbits 
and those with short periods, apsidal precession translates orbital 
differences into angular differences. However, as we discussed in 
Sect. 2.2, the tidal breakup process already spreads the fragments 
along a range of semi-major axes and eccentricities. This means 
that no additional process is required and that the fragments start 
differentially precessing as soon as they form. The inner fragments 
precess more quickly than those on wider orbits, with the timescale 
for complete differential precession of the inner and outer ring (1 
and 2) is given by Eq. 18: 


27 
Ó Tprecess EE EE 
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in the non-dispersive limit where Rast << Rerit For partially un- 
bound orbits, the precession timescale is instead given by 2z/ ogg. 
We visually indicate the differential apsidal precession of an ec- 
centric tidal disc in Fig. 8, where we plot the fragment orbits of a 
100 km asteroid that originates from 3 AU (corresponding to the 
inner boundary of the zone that is not swallowed by the progenitor 
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Figure 8. Differential apsidal precession of fragment orbits from a 100 km 
asteroid originating from 3 AU. The orbits do not cross at t = O (panel a) 
but angular differences increase as a function of time, and at 0.1 Myr (panel 
b) the orbits are clearly seen to cross at both pericentre and apocentre. At 1 
Myr (panel c), both the relative angles and collision velocities have grown 
further and the collision locations are spread out over a wider range in space. 


star during post main-sequence evolution, see Mustill & Villaver 
2012) at 0.1 and 1 Myr. While this is less than the timescale of 
complete orbital precession (~ 20 Myr), the highly eccentric na- 
ture of the tidal disc means that collisions already occur at small 
angular differences. Initially, the orbit crossings are restricted to 
the pericentre and apocentre of the orbit. The relative fragment ve- 
locities scale with the angular differences and increase as a function 
of time. As the angular differences increase, the locations of orbit 
crossings eventually spread out over a wide range in space. 


4.1.1 Differential precession vs coherent precession 


Before we proceed with modeling the collisional grind-down in- 
duced by differential apsidal precession in the next section, we 
should note that not all eccentric astrophysical discs are observed 
to precess at differential rates. For instance, the asymmetric nucleus 
of the Andromeda galaxy (M31) likely consists of an eccentric disc 
(Tremaine 1995; Lauer et al. 2005) that does not show signs of 
differential precession, even though its stellar ages (~ Gyr) far ex- 
ceed the timescale of differential apsidal precession (~ Myr). We 
should, therefore, first investigate whether the same processes that 
act here could cause the orbits of fragments around white dwarfs to 
precess coherently as well. The process that has been suggested to 
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explain the coherent precession of M31 is a dynamical oscillation 
of the eccentricities as a result of self-gravity (Madigan et al. 2018; 
Wernke & Madigan 2019). Physically, any orbits that precess faster 
than the bulk of the disc are pulled back while orbits that lag behind 
are pulled forward. This force torques the orbit and changes its an- 
gular momentum and eccentricity, generating an oscillation. Based 
on the calculations of Madigan et al. (2018), the typical period of 


these oscillations is the secular timescale Le = P (ham), which 


is on the order of 100 orbits for M31. In the case of a disruption 
of an asteroid around a white dwarf, the mass ratio of the central 
object is much larger, however. Even in the case of a disrupting 
Earth-mass planet, the mass ratio (0.6 Mo /Ma) leads to a charac- 
teristic scale around 10° orbits, too long to prevent orbit crossings 
via differential apsidal precession. 

A second option is that differential precession may be inhib- 
ited by continual collisions between fragments. If enough collisions 
occur at small relative angles and velocities, their exchange of an- 
gular momentum could prevent further differential precession sim- 
ilar to what happens in eccentric gaseous discs. For this process to 
be effective, the disc needs to contain a sufficient fraction of small 
particles to generate a large collisional cross section. In the case of 
the tidal discs that form from tidal disruptions around white dwarfs, 
this is unlikely due to the blow-out of the smallest fragments (Sect. 
2.3) and the typically rapid accretion of slightly larger dust grains. 
Hence, it seems safe to assume that differential apsidal precession 
does indeed proceed uninhibited for the larger fragments contained 
in the tidal debris discs around white dwarfs. 


4.2 Gravitational perturbations by a planet 


Besides the gravitational perturbations from the central star, grav- 
itational interactions with surrounding planets can also drive col- 
lisions between the eccentric fragments. Polluted white dwarfs are 
thought to be surrounded by whole planetary systems - including 
potentially gas- and ice giants as well as rocky planets (see Veras 
(2021) for a recent review). Because the tidal disc that forms in 
a disruption event is typically centered around the asteroid's prior 
orbit, interactions with the planet that scattered it will continue to 
perturb the surviving fragments !. 

The effectiveness of this continued scattering likely depends 
on the geometry of the tidal disc. Any fragments whose orbits ei- 
ther cross the path of the planet or whose apocentre is close to the 
planet's semi-major axis will continue to be scattered. The width 
of this chaotic zone (Óacnao, )around the planet has been studied in 
detail in previous works - albeit with far less eccentric orbits - to 


1 
be around 6dchaos = Capi (Œ) * with constant 1.3 <C < 2 (Wis- 


dom 1980; Duncan et al. 1989; Quillen & Faber 2006; Chiang et al. 
2009). If the planet is located at the apocentre of the fragment or- 
bits, this criterion translates to a width of orbit-crossing (Ôa; cross): 
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where we use that the eccentricities of fragments in these tidal 
discs are close to unity. We visualise the range of this scattering 


! During the preparation of this manuscript, we became aware of the con- 
temporary manuscript by Li et al. (2021). These authors have independently 
come to similar conclusions as presented here, based on numerical simula- 
tions of the continued scattering of eccentric fragments. 
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Figure 9. Continued scattering of fragments that cross the chaotic zone of 
a planet that resides at the apocenter of a scattered asteroid. This figure 
depicts the eccentric tidal disc, dividing the fragments between those that 
are likely to be directly re-scattered by the planet (red) and those that now lie 
outside the planet's chaotic zone (green). The top panel (a) shows that most 
fragments from a 50 km asteroid from 3 AU are susceptible to scattering 
(fraction feross). Fragments from a larger 500 km asteroid that originates 
from 30 AU (panel b) are spread bimodally (see Sect. 2) and are mostly 
safe from further scattering. 


for 100 and 500 km asteroids with semi-major axes of 3 and 30 
AU in fig. 9. If the asteroid is small and originates from the inner 
zone, it will form a narrow orbital band and most of its fragments 
lie within the chaotic zone. Larger asteroids or those from further 
out are less likely to be perturbed directly, as the fragment orbits 
begin to follow a bimodal bound-unbound distribution, with up to 
half of the fragments tightly bound to the star with orbits that are 
safe from direct scattering. Nevertheless, slower perturbations of 
these inner orbits are also possible around interior mean-motion 
resonances or via secular perturbations if the planetary orbit is not 
entirely circular, as discussed by Veras et al. (2021). 

Over time, these orbital perturbations will lead to collisions 
between fragments and facilitate faster collisional grind-down. In 
addition, some fragments will be scattered out of the system while 
others are scattered to bound orbits with reduced pericentre dis- 
tances. We suggest that fragment scattering in this way can provide 
a separate avenue for white dwarf pollution. Instead of having to 
lose angular momentum through stellar light, some fragments will 
collide with the white dwarf directly. Others are scattered into the 
sublimation zone, where the vapor pressure of their material be- 
comes significant and they and disintegrate over several orbits. 


4.3 Drag-assisted circularisation via pre-existing material 


Thirdly, fragments can change their orbits by interacting with pre- 
existing gaseous or dusty material inside the Roche radius (Gr- 
ishin & Veras 2019; O'Connor & Lai 2020; Malamud et al. 2021), 
such as has been observed around a substantial minority of sys- 
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tems (Rocchetto et al. 2015; Farihi 2016; Wilson et al. 2019). The 
near-unity eccentricities of fragments mean that their pericentre ve- 
locities are typically on the order of several 100 km/s, yielding 
an extremely violent interaction with any dust or gas that is en- 
countered. In a detailed study, Malamud et al. (2021) showed that 
this interaction can significantly alter fragment orbits. Firstly, the 
fragments lose kinetic energy and become more tightly bound to 
the star. Interestingly, the orbital contraction already becomes sig- 
nificant when the mass of the central compact disc is several or- 
ders of magnitude below that of the tidally disrupted asteroid. The 
smallest fragments circularise the fastest, typically within a few or- 
bits. Larger fragments require additional passages through the disc, 
causing increased orbital differences that again accelerate the onset 
of orbit crossings via differential apsidal precession. In addition, re- 
gardless of whether the pre-existing disc contains dust or gas, each 
fragment passage though the central compact disc erodes away a 
mass similar to the mass that is encountered. Complete circularisa- 
tion of a fragment requires colliding with a similar amount of mass 
as the fragment itself and therefore also leads to its complete disin- 
tegration. In the presence of such a massive central compact disc, 
fragments can accrete onto the white dwarf without the necessity 
for collisional grind-down. 


4.4 Orbital changes due to the Yarkovski effect 


Finally, larger fragments (> 0.1-100 m) may either gain or lose 
angular momentum over time due to the Yarkovski effect (Bottke 
et al. 2006; Veras et al. 2015a,b). The idea is that while a fragment 
orbits the white dwarf, its side that faces the star is more strongly 
irradiated. Subsequent re-emission occurs with a time lag that, cou- 
pled with a rotation, leads to either an accelerating or braking term. 
If the fragment spins sufficiently quickly, its temperature gradient 
smooths out and the term disappears. While the Yarkovski effect 
likely dominates over PR drag in most larger fragments, its derived 
terms are highly dependent on poorly constrained fragment param- 
eters like the spin period (see Veras et al. 2015a, 2019). While a first 
attempt at constraining the spin distribution was made by Mala- 
mud & Perets (2020b), current simulations are not yet able to re- 
solve them for physical fragment sizes. These characteristics of the 
Yarkovski effect currently make a useful inclusion in a collisional 
model unfeasible. However, it is clear that by increasing or decreas- 
ing the orbital momentum of different fragments, the Yarkovski ef- 
fect will both facilitate collisions directly and induce orbital differ- 
ences that accelerate the process of differential apsidal precession. 


5 STAGE IIb: CALCULATION OF ECCENTRIC 
COLLISIONAL GRIND-DOWN 


In this section, we formulate a crude but quantitative calculation of 
collisional grind-down within the highly eccentric tidal discs that 
form after asteroids tidally disrupt around a white dwarf. Our model 
is based on the angular differences that are induced by differential 
apsidal precession (Sect. 4.1), which in turn originate from the or- 
bital spread imparted at the moment of tidal breakup (Sect. 2.2). As 
discussed in the previous section, there are many additional mech- 
anisms that also drive orbital differences between fragments. We 
take differential precession as the sole perturbing process here to 
make the calculation tractable and because it is universally applica- 
ble to tidal discs around white dwarfs, that all start with the required 
orbital spread. 
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Figure 10. Time evolution of the collision rate OM) from our model of 
eccentric collisional grind-down induced by differential apsidal precession 
(see text for details). In the top panel (a), we take asteroids from 10 AU and 
vary their size between 50-500 km. In the bottom panel (b), we take 100 
km asteroids and vary their initial semi-major axis between 3-30 AU. The 
open and filled dots indicate the points where a total of 5096 and 90% of the 
fragment mass has catastrophically collided, respectively. Larger asteroids 
on tighter orbits disrupt to form tidal discs whose fragments collide within 
the shortest period of time. 


5.1 Numerical setup 


We opt for a simple computational approach where we divide the 
fragments into a 2D grid along semi-major axis and fragment size 
and model the collisions with a particle-in-a-box method. The frag- 
ments are assumed to be spread evenly across energy bins, with 
semi-major axes and eccentricities that lie between the bounds 
specified by Eq. 6a-7. Their initial sizes are assumed to follow the 
collisionally evolved distribution of Eq. 11 with œ = 3.5 as a typi- 
cal value, meaning that most of the mass is contained in the larger 
fragments, whereas their surface area is dominated by the smaller 
fragments. As described in Sect. 4.1, the orbits precess at different 
rates depending on their eccentricities and semi-major axes (Eq. 
18), leading to orbit crossings. We estimate the rate of collisions in 
a projected 2D plane, where the collision points can be identified 
at radius reo, and true anomaly 0,4 for any two bins with indices 
i, j from the criterion that rj(co1, i) = rj (col; Øj) = rco, Where the 
velocity v; and distance r; are given by the standard equation of the 
ellipse in polar form: 
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While all collisions exchange some angular momentum, only those 
that are sufficiently violent lead to the (partial) destruction of frag- 
ments. Accounting for the angular momentum changes in sub- 
catastrophic collisions is not possible with our method due to the 
large number of additional spacial bins it would create. Hence, we 
do not include these less violent collisions in our calculation and 
only focus on the collisional grind-down from catastrophic colli- 
sions. The required specific energy for catastrophic fragmentation 
Q-l(ü- Dik m is known as the dispersion threshold and can be 
estimated as a scaling relation (Durda et al. 1998; Benz & Asphaug 
1999): 


O* = Q,R " + QR”, (22) 


where the first term accounts for the material strength of a fragment 
and the second term corresponds to the gravitational binding energy 
that has to be overcome. We take the constants a = 0.3, b = 1.5 
and Q, = 6.2 - 10" erg/g, Qs = 5.6 - 10^? erg/g from Löhne et al. 
(2008) and Wyatt et al. (2011). In our simulations, fragments of 
different sizes that reside in the same orbital bin experience the 
same collisional velocities when they cross paths with fragments 
in other bins. This means that Eq. 22 also directly specifies the 
minimum fragment size Rj crit that can catastrophically collide with 
a fragment of size R;, depending on its collisional velocity: 
1 


agis 
Rj crit = E3 Ri. (23) 


col 


Generally, only fragment pairs with size ratio’s below two orders 
of magnitude are found to collide catastrophically. 

The collision rate Pj of two bins with fragment sizes R;, Rj and 
Nj, Nj fragments can be estimated from a standard particle-in-a-box 
approach as: 


X Thox,i > < Toox > NiNjOijVcol,ij 


Pj= , (24) 


Torb,i Torb.j Vox 

where oi; = z(Ri TR is the collision cross section, Vpox is the 
volume of the collision box and < Tyox,i 7, € Tbox,j > are the av- 
erage periods of time that fragments of bins i and j spend there. 
We approximate these collision volumes as locally straight boxes 
with volume Viox ~ ll; Hsin(05;) where ou = |vi x vi|/(|vi||vi|) is 
the angle between the orbits and H is the local height of the tidal 
disc. We estimate this height based on a typical orbital inclination 
i as H = 2irgo, = 2Rastřcol/rg (see Sect. 2.2), assuming that the 
inclination remains constant over time. The widths J; and J; cancel 
from the time spent in the box «Tio, ; >= li/|Vi|, which gives: 


| ZNN (Ri +R;)? — vi 
Toro iTorb j H | Vi x il 


d (25) 


The total collision rate of bin i is set by the binned sum over P; j: 


Jerit,max 
P= } Pj, (26) 

Jcrit,min 
which we evaluate numerically. Because we only track catastrophic 
collisions in this scheme, we register the first catastrophic en- 
counter for any fragment and remove the pair involved in the col- 
lision for the rest of the simulation. Although this approach does 
not incorporate grind-down that can result from second genera- 
tion fragments, larger fragments generally require more time to 
catastrophically collide so our approach can still provide order-of- 
magnitude estimates. More importantly, even this crude model can 
illuminate important trends and biases that are inherent to the col- 
lisional phase, which we will describe in the next subsections. 
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Figure 11. Timescales for the collisional grind-down of fragments in eccen- 
tric tidal discs (t4), calculated with our numerical model based on differ- 
ential apsidal precession (see text for details). The grind-down timescale is 
defined as the time required for half of the mass to catastrophically collide. 
Fragments within tidal discs that form from large asteroid progenitors on 
tight orbits grind down within the shortest period of time. 


5.2 General trends in eccentric grind-down rates 


We explore the most important trends of eccentric fragment colli- 
sions in Fig. 10 where we simulate the grind-down of fragments 
from a range of asteroid progenitors. In the top panel (a), we take 
asteroids with fixed semi-major axes of 10 AU but with sizes be- 
tween 50-500 km. In the bottom panel (b), we instead take fixed 
asteroid sizes of 100 km but vary their semi-major axis between 
3-30 AU. The first thing to note is that their rates of grind-down 
follow similar temporal shapes. Initially, no fragments cross paths 
and the rate collisions begins at zero. As the fragment orbits con- 
tinue to precess at different rates, some orbits begin to cross but the 
relative velocities are low and only similar-sized fragments catas- 
trophically collide. After 10? — 10^ years, these relative velocities 
have increased sufficiently that catastrophic collisions occur for a 
wide range of relative fragment sizes and the collision rate shoots 
up. When around half of the fragment mass has collided (open 
dots), the rate decreases from its peak value and continues to drop 
as fewer and fewer intact fragments remain. 

The most important variable in the rate of grind-down is the 
size of the asteroid progenitor. Larger asteroids provide more frag- 
ments that can collide (factor of Ria), which additionally means 
that each fragment has more other fragments that it can collide 
with (also a factor of R3„). Together, this predicts a scaling of the 
grind-down rate of Moo, œ RÊ, similar to what we find in our sim- 
ulations. In Fig. 10a, we find that the peak accretion rate increases 
from ~ 108 g/s for a 50 km asteroid from 10 AU to ~ 10? g/s for a 
500 km asteroid with the same semi-major axis. This slightly flatter 
scaling (~ RÀ, instead of œ RÓ.) is due to the more disperse disc 
that forms when a larger asteroid disrupts (see Figs. 3 and 4). Some 
fragments of the 500 km progenitor are placed on unbound orbits 
and are ejected form the system, while others are just placed on 
wider orbits that take longer to collide. This general inverse scaling 
of the rate of grind-down with the semi-major axis of the asteroid 
is also shown in Fig. 10b, where a 100 km asteroid from 3 AU is 
found to have a peak collision rate that exceeds the collision rate of 
an asteroid from 30 AU by around an order of magnitude. 
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Figure 12. Cumulative time where the rate of dust production from colli- 
sional grind-down lies above a certain rate t(Mzo, > Mx). The values corre- 
spond to calculations with our collisional model (see text for details). Larger 
asteroids produce more fragments, leading to more rapid grind-down and 
higher peak collision rates. Asteroids that originate from the inner zone of 
planetary systems (solid lines, 3 AU) generate higher peak collision rates 
compared to those on wider orbits (dashed lines, 10 AU) but they spend less 
time accreting at more moderate rates, making their pollution less likely to 
be observed if the subsequent accretion of small dust is sufficiently rapid. 


5.3 Event lifetimes and peak accretion rates 


In order to compare the results of our simulations to observationally 
inferred accretion event lifetimes, we simulate a grid of 400 tidal 
discs corresponding to asteroid sizes between 50-500 km that orig- 
inate between 3-30 AU. Our results are shown in Fig. 11, where we 
indicate the time required to grind down half of the fragment mass 
(labeled 1,51). As explained in the previous subsection, our model 
yields a large variation in collisional timescales based on the aster- 
oid size and its semi-major axis. Rather than try to produce an exact 
fit to the DAZ population from a crude model, we investigate the 
observability trends based on these parameters. The shortest grind- 
down timescales that we find are on the order of 10? years for as- 
teroids that are several hundreds of kilometers in size. When the 
asteroid size is reduced to 50 km, the timescale increases by orders 
of magnitude to 10° — 10" yr, depending on the orbital separation 
of the asteroid. From an observational perspective, typical accretion 
lifetimes of material around white dwarfs can most directly be in- 
ferred from differences in detection rates of DAZ and DBZ stars. In 
the pioneering study by Girven et al. (2012), the typical timescales 
are inferred from the difference between DAZ and DBZ pollution 
rates to be between 10^ — 10° yr. More recently, this analysis was 
revisited by Cunningham et al. (2021) with updated photospheric 
modeling, which yielded an order of magnitude longer timescales 
between 10? — 107 yr. A different approach was taken by Harri- 
son et al. (2021), who estimated typical accretion event lifetimes 
around 107 yr based on a Bayesian analysis of the photospheric 
composition. When compared to our simulation results, lifetimes 
around 104 — 106 yr match our computed grind-down timescales 
for 100-400 km asteroids, whereas we only find longer lifetimes of 
107 yr when the asteroids are smaller than 50 km. 

The highest inferred on-going accretion rate to date is 10? 
g/s for the most metal-rich DAZ (Farihi et al. 2012). If convective 
overshoot is accounted for, the inferred rates could increase by an- 
other order of magnitude (Cunningham et al. 2019). In Fig. 12, we 
examine whether such high accretion rates can be generated by the 
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grind-down of asteroids with sizes between 100-300 km that origi- 
nate from either 3 or 10 AU. The figure shows that this is indeed the 
case when any produced dust is rapidly accreted, as all lines show 
at least some period of grind-down rate beyond 10?? g/s. Within 
the context of our model, grind-down rates beyond 10!! g/s require 
asteroids with sizes that exceed 200 km. 

The most extreme outliers of the inferred accretion rates have 
only been observed for DBZ stars. Farihi et al. (2012) identified six 
objects with inferred accretion rates above 1010 g/s, with the high- 
est observation implying a record rate of 10!!? g/s. Due to their 
longer sinking timescales, these are not observations of on-going 
accretion rates but represent an average over a longer time period. 
Based on their findings, Farihi et al. (2012) suggest that short peri- 
ods of violent accretion that take between 10-100 years must occur 
at times to explain the difference between DAZ and DBZ upper 
values. Within the context of our model, the dichotomy between 
the DAZ and DBZ upper values arises naturally. More massive as- 
teroids produce more fragments that each collide within a shorter 
period of time, causing the accretion of the biggest asteroids to oc- 
cur in short and intense bursts, provided that the accretion of the 
produced dust is sufficiently fast. In contrast, the grind-down of 
fragments that originate from smaller asteroids occurs over longer 
periods of time, making them more commonly detected in DAZ 
pollution, which measures on-going accretion. 


5.4 Monte-Carlo analysis of accretion rate distribution 


If the dust that is produced in the collisional grind-down of larger 
fragments is accreted quickly, the collision rate directly translates 
into an accretion rate onto the star. In this subsection, we explore 
the trends and biases that this would imply for observations of pol- 
luted white dwarfs. We perform Monte-Carlo simulations for a few 
different asteroid populations, where we scatter and tidally disrupt 
a total mass of 2.4 -10?4 g in asteroids over a period of 1 Gyr, which 
is the equivalence of the mass contained in the Solar System’s main 
belt (Pitjeva & Pitjev 2018). While this amount may seem exces- 
sive, it is also equivalent to only 296 of the Kuiper belt (Fraser et al. 
2014) or 1% of the mass of the debris disc around € Leporis (Chen 
& Jura 2001). Spread over 1 Gyr, it amounts to average accretion 
rate between 3.9 — 7.8 - 107 g/s. 

We contrast the statistical distribution of our calculated grind- 
down rates with the DAZ sample of Koester & Wilken (2006) and 
Zuckerman et al. (2003), which was compiled and used for the 
same purpose by Wyatt et al. (2014). This sample is specifically 
selected for this purpose because the stars it contains are randomly 
chosen based on being nearby and bright, and they are not biased 
in terms of the presence or absence of metals. The sample contains 
a total of 534 DAZ systems from a Keck and SPY survey, with 38 
CA detections and 298 upper limits on CA. We only use the 467 
systems that have estimated sinking times less than 104 yr (from 
the estimate of Koester 2009), so that the sample reflects on-going 
accretion rather than an average over past accretion events. In or- 
der to make a uniform comparison across the sample, we follow 
Wyatt et al. (2014) and calculate the total accretion rate based on 
the Ca accretion rate with a scaling that matches bulk Earth, noting 
that true accretion rates may vary from this assumption, particu- 
larly in Ca-rich bodies, by orders of magnitude. Using the sinking 
times from Koester (2009), the measurements imply accretion rates 
between 109-! — 10? g/s. Identical to the methodology of Wyatt 
et al. (2014), the cumulative probabilities are calculated from the 
sub-sample of systems Ns (Mx) where accretion could have been 
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Figure 13. Cumulative probability distribution of grind-down rates (colored 
curves) compared to an unbiased young DAZ sample of inferred accretion 
rates (grey, compiled by Wyatt et al. 2014), plotted with a 1 o uncertainty 
band. The grind-down curves correspond to Monte-Carlo simulations of our 
eccentric collision model (see text for details). In this plot, we examine sys- 
tems where the equivalence of the main belt (1.2 — 2.4- 1074 g) is scattered 
in over a period of 1 Gyr. The colored lines correspond to mono-size as- 
teroid distributions of 100 km (blue), 200 km (green) or 300 km (red). The 
asteroids have semi-major axes of 3 AU (solid) or 10 AU (dashed). 


detected at that level. The 1 o errors are calculated from binomial 
statistics (Gehrels 1986). 

We plot the DAZ sample together with the results of our 
Monte-Carlo simulations in Fig. 13. The model presented here pre- 
dicts accretion rates based on two parameters, namely the asteroid 
size and its initial semi-major axis. As explained in the previous 
subsection, larger asteroids produce more fragments that each col- 
lide within shorter periods of time, and hence appear as short but in- 
tense bursts of accretion. This can be compared to the slower, stead- 
ler accretion from smaller asteroids (see Fig. 10). When the same 
total mass in asteroids is scattered towards the star, a population 
containing larger asteroids leads to a steep cumulative distribution 
of the accretion rate, with few white dwarfs found accreting at the 
highest rates, whilst a population of smaller asteroids leads to a flat- 
ter distribution with more white dwarfs found at typical accretion 
rates but lower peak values. These trends are visualised in Fig. 13 
from the comparison between curves with the same line style (semi- 
major axis) but different colors (size). Each of the curves within the 
same panel represents the same total mass in asteroids scattered in, 
with the red curves representing accretion of the biggest asteroids 
(300 km) and the blue corresponding to the smallest asteroids (100 
km). While the bigger asteroids have higher peak accretion rates 
(see also Fig. 12), their cumulative probabilities of accretion at de- 
tectable levels (> 107 g/s) are as much as an order of magnitude 
lower. 

The observability bias based on semi-major axis proceeds 
along the same lines. Asteroids with wider orbits prior to their tidal 
disruption release their fragments along more spread-out discs. As 
a result, their fragments each take longer to complete an orbit and 
require more time to catastrophically collide. This means that the 
fragments from wide asteroid progenitors spend more time produc- 
ing dust at lower rates, at the cost of having a reduced peak value 
(see Fig. 12). In the Monte-Carlo simulations of Fig. 13, we visu- 
alise the difference with solid and the dashed curves of the same 
color, which correspond to the accretion of identical asteroids with 
semi-major axes of 3 AU and 10 AU, respectively. The effect of or- 


bital separation exposes a second bias towards detecting on-going 
accretion at moderate values ( 10? g/s) by asteroids that originate 
from wider orbits. 


6 STAGE III: INFRARED EMISSION FROM 
ACCRETING FRAGMENTS OR DUST 


A substantial minority of accreting white dwarfs exhibit detectable 
infrared excesses (Rocchetto et al. 2015; Farihi 2016; Wilson et al. 
2019). While the simplest canonical model computes emission 
from opaque dust on circular orbits (Jura 2003), elliptical discs 
around white dwarfs have been modeled in response to obser- 
vational clues to account for the reduced infrared excess around 
young stars (Dennihy et al. 2017) and also to explain the variability 
of the emission Nixon et al. (2020). The elliptical nature of the ma- 
terial is furthermore suggested by Doppler tomography of gas near 
the Roche radius (Manser et al. 2016; Steele et al. 2021). From 
a physical perspective, it is similarly clear that the fragments are 
indeed expected to be released on orbits that are initially highly ec- 
centric (see Sect. 2.2) and then slowly circularise over time. In our 
modeling, we compute the emission that is generated during the 
accretion of dust from such a highly eccentric structure. We do not 
follow the typical assumption that the disc is opaque to stellar light. 
Instead, we compute the distribution of dust and its re-emission un- 
der the a-priori assumption of an optically thin environment, the 
validity of which we then discuss after we calculate the radial opti- 
cal depth. 

With our model, we show that that typically inferred accretion 
rates can indeed produce detectable infrared excesses when grains 
circularise in an optically thin disc via PR drag alone, but that the 
disc only remains optically thin if the orbital inclinations of the 
dust are increased over time. We then evaluate the results within 
our three-stage model and discuss how the disc geometry and cir- 
cularisation speed can effect the produced infrared excesses. 


6.1 Infrared emission from a collision-less tidal disc 


Before we calculate the infrared emission by dust grains from col- 
lisional grind-down, we first evaluate how much emission is pro- 
duced just after a tidal disc forms and how this changes as frag- 
ments circularise. The emission produced by newly formed tidal 
discs was suggested by Nixon et al. (2020) as the source of ob- 
served infrared excesses based on normalised emission profiles. 

As explained in Sect. 2.2, the fragments are initially placed on 
a range of orbits depending on the size of the asteroid, that share 
their pericentre but differ in their apocentre (see Fig. 4). In this cal- 
culation, we simplify the situation somewhat and assume that the 
fragments initially occupy a single orbital ring (ao, eo). given by the 
original orbit of the asteroid. Along this ring, the fragments experi- 
ence different temperatures, which we approximate with radiative 
equilibrium: 


Nie 


T(r) =Two ( - ) (27) 
Rwp 

For simplicity, we ignore the super-heating of the smallest parti- 

cles in this model (Chiang & Goldreich 1997; Rafikov & Garmilla 

2012). The Planck emission per unit mass (.Zy) depends on the size 

of the fragments, which determines their ratio of area (A) to mass 
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Figure 14. Simulated emission per unit surface area at 10 um from fragments accreting onto G29-38 (0.59 Mo, 11240 K, 17.5 pc) via collision-less PR drag, 
computed assuming optically thin properties (See Sect. 3.1 and Sect. 6.2 for details on PR drag and this emission). The three panels correspond to snapshots at 
different times after the tidal disruption event of an asteroid from 3 AU and the colors are normalised individually per panel. Just after the tidal disruption event 
(a), the fragments are still highly eccentric and their emission is minimal. The smallest fragments accrete just after 1900 yr (b), when their more contracted 
orbits begin to generate far more emission. The final panel corresponds to 10? yr, when all fragments smaller than 100 um have already reached the star. The 
disc contains an elliptical inner gap because fragments are not yet fully circularised when they enter the sublimation distance (white dashed circle). The Roche 


radius is indicated with a red dotted circle, just outside the sublimation zone. 


Rast = 750 km 
10! ji 
> Q0 
E 10 
u? 
10-71} SO S 
--- t=1900 yr WK 
10-2 
1071 109 10! 102 
À [um] 


Figure 15. Simulated emission spectra from a 750 km asteroid disruption 
around G29-38 computed by our model (see text), followed by collision- 
less and optically thin orbital contraction via PR drag. The dotted line in- 
dicates the contribution from the star, with the dash-dotted (t — 0), dashed 
(t — 1900 yr) and solid (t — 10? yr) lines including the emission from the 
fragments at different times. The points indicate the emission from the sys- 
tem as observed in several surveys (largely compiled by Farihi et al. (2014), 
original references in the text). 


(M): 
Fy(R,7) = 250) (28a) 
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where By is the Planck function at a given frequency v. In this sim- 
plified collision-less test case, fragments evolve towards the star 
via PR drag as a function of time, with their semi-major axes and 
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eccentricities slinking at rates given by Eqs. 13a and 13b (Veras 
et al. 2015a,b). At any time f, fragments with different sizes R 
have already shifted away from their initial orbit and are now lo- 
cated at their own orbits (a(t),e(t)), which are most contracted 
for the smallest fragments. We take a linear grid of angular bins 
with size dO between 0 and 27, where fragments spend a time 
dt = d0/6(@). The average emission from a particular orbit per 
unit mass (Fy orbit) is then weighed by the time spent at given an- 
gle 0 as: 


" 1 Z dë 
Fy osi (R, t) = Des Y Ful Rit, er (29) 
1) @=0 


where P is the fragment's orbital period. The total emission (Fy pr) 
from the asteroid fragments is given by the sum over all the rings 
occupied by the different fragment sizes at time f: 


Rmax dM 
Fypr() = Y, Fv ovi (Rt) Ae (PAR, (30) 
Rmin 


where dM /dR is defined by the fragment size distribution. As ex- 
plained in Sect. 2.3, we take a truncated power law with the default 
exponent @ = 3.5. The definite lower limit is set by the blow-out 
size (Eq. 10) and the upper limit is set at 1 km, corresponding to a 
rough indication of the breakup limit from the tidal force (Eq. 4). 
We visualise the results by calculating the expected emission 
around G29-38, a 11240 K DAZ white dwarf at 17.5 pc with a 
well-documented infrared excess (Xu et al. 2018), first identified 
by Zuckerman & Becklin (1987). We compare our results to multi- 
wavelength fluxes, for which we take a supplemented version of 
the sample compiled by Farihi et al. (2014) with observations from 
Tokunaga et al. (1990), Reach et al. (2005), and Farihi et al. (2008). 
For our comparison here, we disrupt a large, 750 km asteroid with 
initial semi-major axis at 3 AU and let the fragments circularise in 
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a collision-less, optically thin manner via PR drag. Without colli- 
sions, this scenario produces a peak accretion rate just above 10? 
g/s that occurs around 1900 yr after the disruption event. We show 
the geometry of the emitting fragments around the white dwarf in 
Fig. 14, plotted at three different times at an emission wavelength 
of 10 um. Panel a (t = 0) corresponds to the moment after dis- 
ruption, assuming that the fragments are spread uniformly over the 
orbit, which is known to only require a few orbits (Veras et al. 2014; 
Malamud & Perets 2020a; Li et al. 2021). Panel b is plotted at 1900 
yr, just before the first fragments sublimate and accrete, which hap- 
pens when their pericentre crosses the sublimation distance (white 
dashed circle), assumed to be located at 1500 K, corresponding to 
the temperature where the vapor pressure of silicates becomes sig- 
nificant and they begin to sublimate in vacuum (van Lieshout et al. 
2014). Panel (c) corresponds to 10? yrs after the disruption event, 
when fragments smaller than 100 um have already accreted. In all 
three cases, the non-circular nature of the material is evident. Most 
emission at 10 um occurs near the orbit's pericentre, or just away 
from it in the fragment's most eccentric initial state. Because the 
fragments enter the star's sublimation radius before they are en- 
tirely circular, the central gap has the form of an ellipse rather than 
a circle. The eccentricity of this ellipse increases for higher stellar 
temperatures and for lower sublimation thresholds. 

In Fig. 15, we plot the emission spectrum of the fragments at 
the same snapshots in time. The first thing to note is that the emis- 
sion varies greatly with time. The initial fragment orbits have semi- 
major axes similar to that of the original asteroid (3 AU here), and, 
therefore, have long orbital periods with most of their time spent at 
a distant apocentre. As a result, the emission from the orbit's peri- 
centre is initially minimal, with emission only picking up at longer 
(~ 100 um) wavelengths. Based on these calculations, we argue 
that newly formed tidal discs provide insufficient infrared emission 
to explain the observed excesses or variation, as was suggested 
by Nixon et al. (2020). We find that the emission from the disc 
does rapidly increases as the fragments circularise, peaking in this 
collision-less test case at the moment that the smallest fragments 
begin to reach the star. The emission then gradually decreases again 
as the smallest fragments reach the white dwarf, with the emission 
spectrum maintaining a similar shape. During the circularisation 
process, the effective temperature of the emission slightly increases 
over time as fragments contract their orbits and spend more time at 
strongly irradiated locations. 


6.2 Emission from collisional dust production 


We can slightly modify the calculation of the previous subsection to 
calculate the emission produced by the accretion of small dust that 
is produced via collisional grind-down of larger fragments. With 
the same assumption of optically thin contraction via PR drag, we 
can calculate the average emission per unit mass (.Zy) of a frag- 
ment from the moment of its production to its accretion by averag- 
ing over the different orbits (a(t), e(t)) that a dust grain occupies on 
its trajectory towards the star: 


Fy (R) = Fy orit(R,t’ dt’ (R,t'). (31) 


The time intervals dt’ during which the different orbital bins are oc- 
cupied follow from Eqs. 13a, 13b. We again sum the contributions 
of different fragments in the size distribution to obtain the total 
emission (Fy pg). but now for a given accretion rate rather than at a 
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Figure 16. Simulated emission per unit surface area at 10 um from constant 
dust production and accretion at a rate of 107 g/s via PR drag, computed 
assuming optically thin properties. The top panel (a) shows the resulting 
emission without apsidal precession and the bottom panel (b) shows the 
axisymmetrically averaged emission. The disc contains an elliptical inner 
gap because fragments are not yet fully circularised when they enter the 
sublimation zone (white dashed circle). The Roche radius is indicated with 
a red dotted circle, just outside the sublimation zone. When axisymmetri- 
cally averaged, the structure becomes visible as a ring-like structure mostly 
contained within the Roche radius. 


particular time: 


2o nn ER dM 
Fy pu (M) = ), y (Rytacc(R) Fe (R)AR, (32) 
Rmin 


which, for a given accretion rate, notably yields a result that is inde- 
pendent of the fragment size distribution if the fragments are mod- 
eled as black bodies in thermal equilibrium ?. We show the geome- 
try of the emission in Fig. 16, which, as expected, shows a similarly 
asymmetric shape as the time evolution discussed in the previous 
subsection, with most of its emission originating from a narrow 
band near the pericentre of the orbit. If the collisions are induced 


2 The PR accretion time scales as face e R, while the area-to-mass ratio 
scales as A/M œ R-, These contributions cancel out, meaning that the 
emission in a steady-state of dust production and accretion by PR drag is 
independent of R if the particles are modeled as perfect absorbers/emitters. 
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Figure 17. Simulated emission spectra of G29-38, with dust accretion via 
PR drag at three different rates. The dotted line indicates the contribution 
from the star, with the dash-dotted (107 g/s), dashed ( 108 g/s) and solid 
(10? g/s) lines including the emission from circularising dust at different 
accretion rates. The points indicate the observed emission from the system 
in several surveys (largely compiled by Farihi et al. (2014), original refer- 
ences in the text). 


by the differential apsidal precession of fragment orbits, the ac- 
tual geometric shape of the disc will be different. We approximate 
this state in panel (b) of the same figure, where we neglect poten- 
tial collisional circularisation from inelastic collisions, yielding an 
axisymmetrically averaged structure when the fragment orbits are 
completely differentially precessed. The pericentre now appears as 
a bright narrow band up to the Roche radius, surrounded by a less 
luminous zone. 

We plot the emission spectrum, which is not affected by asym- 
metric averaging, in Fig. 17, corresponding to accretion onto G29- 
38 at three different rates. Because the emission is proportional to 
the assumed accretion rate, the implication of this model is that an 
excess in the system’s infrared emission is visible whenever the PR 
accretion rate is sufficiently high. While such an excess has indeed 
been observed for some rapidly accreting systems, the majority of 
polluted white dwarfs do not follow this trend. Hence, we must con- 
clude that even if PR drag can supply high accretion rates of small 
dust, the majority of systems likely accrete their material in a dif- 
ferent manner. In the next subsection, we will evaluate the optical 
depth of the accretion via PR drag, which provides an additional 
constraint to the dust accretion process via discs with small incli- 
nations. 


6.3 Radial optical depth of the tidal debris disc 


In the previous calculations regarding accretion via PR drag, we as- 
sumed that the dust and fragments around the white dwarf form an 
optically thin disc. Here, we investigate the validity of this assump- 
tion and compute the optical depth in the radial direction, where the 
disc is most opaque. We follow a similar procedure as in the pre- 
vious subsection and begin by computing the radial optical depth 
contribution per unit mass (.7) as a function of the true anomaly 
for a given orbit: 


mR? dt 


FT (R,t,0) = Hg P (33a) 


3 
Apgas RH (R,t, 0)P(R,t)v(R,t, 0)’ 


(33b) 
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where P(R,t) is the period of a fragment at its given orbit. For this 
calculation, we assume the same constant inclination imparted at 
the moment of breakup of i = Rast/rg, which is around 107^ for a 
100 km asteroid (see Sect. 2.2). To obtain the optical depth contri- 
bution per unit mass (.7 (LR. @)) for a fragment size bin, we again 
sum over the different orbital rings that are occupied during the 
orbital contraction, evaluated separately for every angle: 


facc (R) 
—— Z (R,t',0)dt'(R,t'). (34) 
e à, 705a) 
The total optical depth is then determined by the sum over the size 
distribution: 


Z(R,0)— 


. Rmax Zeg dM 
Tr(M,0) — Y T (R,8)tcc(R) e (R)dR. (35) 
Rmin 


We are primarily interested in estimating the optical depth at peri- 
centre because that is where the fragments are exposed to the 
most stellar light, even when corrected for travel time (Veras et al. 
2015a,b). When we perform the calculation with the parameters of 
G29-38 for the star and assume a 100 km asteroid, we find that the 
disc becomes radially opaque at its pericentre when the accretion 
rate is increased beyond 107 g/s. This result is not altered much 
when the material in the disc is axisymmetrically averaged, mean- 
ing that the optical depth at pericentre is comparable to the average 
across the disc. 

While stellar rays from more vertical directions can still reach 
the dust when it is radially opaque, our calculation suggests that 
the assumption of perfect optically thin accretion begins to fail at 
higher accretion rates. If the inclination of the dust remains small, 
the amount of stellar light that can reach the dust grains declines 
and circularisation rates slow down (see also the models by Rafikov 
2011b,a). However, if the inclination of the tidal disc does increase 
from its initial value, for instance due to interactions with neigh- 
boring planets (Li et al. 2021), the disc can remain entirely opti- 
cally thin and high dust accretion rates by PR drag are possible. 
This scenario could explain the few cases where infrared excesses 
are observed. Finally, several white dwarfs show the remarkable 
combination of both rapid ongoing accretion and no infrared ex- 
cess. In the scenario of grind-down and dust accretion, we find that 
this requires faster dust circularisation and accretion than is possi- 
ble by PR drag alone, as was also suggested by Bonsor & Wyatt 
(2010). This finding points to the importance of additional circular- 
ision processes that could involve gas drag (Malamud et al. 2021) 
or the recently suggested mechanism of Alfvén-wave drag (Zhang 
et al. 2021). 


7 DISCUSSION 


In this work, we studied how material accretes onto white dwarfs 
from their surrounding planetary systems and how this relates to 
observational quantities. Our baseline scenario begins with the tidal 
disruption of an asteroid close to the white dwarf, which forms 
a highly eccentric tidal debris disc. The fragment orbits are then 
perturbed via various processes, including differential apsidal pre- 
cession, causing the larger fragments to collide on their eccentric 
orbits until only dust remains. In the final stage, the dust accretes 
onto the star by drag forces. Our suggested scenario can produce 
accretion rates as high as those observed (> 10!! g/s) from the dis- 
ruption of > 200 km asteroids. Both the presence and the absence 
of infrared emission can be explained depending on the rate of dust 
in-spiral and accretion, with drag rates faster than PR-drag, such 


16 Brouwers Et Al. 


as via additional gas or Alfvén-wave drag required to explain the 
absence of infrared emission in the majority of white dwarfs. 

Nevertheless, we propose that no single accretion scenario ex- 
plains the pollution of all white dwarfs. A range of different ac- 
cretion channels are likely applicable depending on the properties 
of both the central star and the accreted object. In this discussion, 
we will present a selection of different routes to white dwarf pollu- 
tion, which are summarised in the form of a road-map in Fig. 18. 
We visualise this road-map as a series of forks that split into differ- 
ent accretion channels according to specified physical criteria. We 
also discuss the observational characteristics that belong to each of 
these channels wherever they are sufficiently well understood, not- 
ing that further detailed modelling is required in many cases. After 
presenting our road-map, we will discuss the limitations of our cal- 
culations relating to our suggested main path and suggest areas for 
improvement in our model. 


7.1 Fork (1): Direct asteroid impact, ejection or tidal debris 
disc formation 


The first step towards white dwarf pollution starts with mass loss 
from the central star, which widens the chaotic zone around planets 
(Bonsor et al. 2011; Mustill et al. 2018) and can destabilise tightly- 
packed planetary systems (Debes & Sigurdsson 2002; Maldonado 
et al. 2020, 2021). Nearby asteroids become subject to scattering 
from close encounters or strong perturbations in mean motion res- 
onances. In principle, scattering events can have four possible out- 
comes, which we show at the first fork (1) of Fig. 18. Most com- 
monly, the asteroid’s eccentricity or semi-major axis are only al- 
tered slightly and the asteroid continues on its way until it is scat- 
tered again. In a chain of scattering events, the asteroid can eventu- 
ally attain such a high eccentricity that it enters the Roche radius of 
the star and it disrupts into an eccentric tidal disc (red arrow). This 
corresponds to our suggested baseline model. Alternatively, the as- 
teroid could either be scattered outwards and enter the influence of 
outer planets, become completely unbound from the system, or di- 
rectly hit the surface of the white dwarf if its pericenter distance 
becomes sufficiently small. 

This last possibility of a direct asteroid impact is worth men- 
tioning as a separate channel of accretion, studied in detail by 
Brown et al. (2017) and McDonald & Veras (2021). It is the sim- 
plest method of mass accretion, as material almost instantly enters 
the star’s photosphere. This near-instant accretion prevents any de- 
tectable infrared excess and also restricts the detection window of 
the pollution itself to a few sinking timescales of metals in the pho- 
tosphere, making direct asteroid impacts unlikely to be detected 
in young DAZ stars. In any case, direct impacts should be rare 
events. Not only is the Roche radius significantly larger than the 
white dwarf itself, but Veras et al. (2021) show that most low mass 
(terrestrial) planets provide small eccentricity kicks that marginally 
push asteroids into the white dwarf’s tidal disruption zone. 


7.2 Fork (2): Sublimation, continued scattering by a planet 
or orbital perturbations and collisions 


We continue along the main channel of our road-map (visualised 
as red in Fig. 18) with the formation of an eccentric tidal disc. Fol- 
lowing the asteroid’s disruption, its fragments can evolve in a num- 
ber of ways as indicated at the second fork (2). If the temperature 
at the disc’s pericenter exceeds the sublimation threshold, its frag- 
ments quickly turn into gas. This is always the expected outcome 
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Figure 18. Road-map outlining potential routes for planetary material to 
arrive in the atmospheres of white dwarfs. Our suggested main route (red 
arrows and boxes) begins with the injection of an asteroid into the stel- 
lar Roche radius, followed by a tidal disruption event, orbital perturba- 
tions, collisional grind-down, and finally dust accretion. Alternative accre- 
tion channels are shown in purple with physical selection criteria at five 
numbered points. The detectable characteristics of these different accre- 
tion channels are listed in green, provided that they are sufficiently well- 
constrained. 
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around hot white dwarfs (= 20.000 K), which are sufficiently lumi- 
nous that rocky material begins sublimating at distances outside the 
star’s Roche radius (Bonsor et al. 2017). However, in rarer cases, 
fragments can also sublimate around less luminous stars if their 
pericenter lies deep in the Roche sphere or when the fragments are 
composed of volatile components, which always sublimate within 
the Roche radius (Steckloff et al. 2021). Although it is clear that 
fragment sublimation leads to a distinct channel of accretion, fur- 
ther work is required to determine the full details. If the subse- 
quent viscous evolution of the gas is sufficiently rapid, the gas can 
quickly accrete onto the star after it is produced, leading to a sce- 
nario of pure gas accretion with potentially detectable gas emission 
lines but no infrared emission. If the gas does not circularise suffi- 
ciently within a single orbit, it likely re-condenses on its way back 
towards apocentre, in which case dust exterior to the Roche limit 
could produce detectable infrared emission. In any case, this sce- 
nario is likely characterised by a high peak accretion rate as small 
fragments quickly sublimate and the presence of gas only adds to 
the circularisation and accretion speeds of remaining solids. 

For those tidal discs where sublimation does not occur, even 
at pericentre, we consider two further possibilities. If the asteroid 
was scattered by a planet, this planet could potentially re-scatter 
the disrupted material. This would occur for those fragments whose 
apocentre continues to approach that of the planet (See sect. 4.2). 
Otherwise, the fragments will evolve according to further orbital 
perturbations, including apsidal precession, leading to collisions. 
We separate these two scenarios because they potentially lead to 
different observational outcomes, as discussed next. 


7.3 Fork (3): Outcome of continued scattering by a planet 


We visualise the possible outcomes of continued scattering after 
fork (3) of Fig. 18. Planets can either scatter fragments outwards, 
such that they are ejected, inwards, such that they graze the star 
with reduced pericenter, or just perturb their orbits, leading to fur- 
ther collisional evolution. The likely outcome depends in part on 
the size of the asteroid and on its semi-major axis. Since larger as- 
teroids and those that originate from an outer belt disrupt into wider 
tidal discs, strong scattering by a planet is unlikely for most of 
their fragments, such that we envisage continued collisional grind- 
down as the most likely pathway in these cases. This preference 
towards collisional evolution is further amplified for large asteroids 
due their fragments shorter collisional time-scales. For smaller as- 
teroids, collisions require more time and many fragments remain 
on planet-crossing orbits, making them instead susceptible to con- 
tinued strong scattering. We predict that as the sublimation zone is 
significantly larger than the white dwarf, and the bodies are deep 
in the white dwarf’s potential, inward scattering of fragments typ- 
ically leads to their sublimation rather than a direct impact. In this 
sense, the accretion channel via planetary scattering might proceed 
similarly to the sublimation-based channel mentioned earlier. Un- 
derstanding the full details will require further work, where a de- 
tailed understanding of gaseous evolution and condensation on the 
highly eccentric orbits will be crucial. 


7.4 Fork (4): Rapid circularisation or collisional grind-down 


We continue our suggested main road at fork (4) of Fig. 18 with the 
evolution of fragments that are not sublimated form the heat of the 
central star nor scattered by a planet. These fragments nevertheless 
have their orbits perturbed via various processes (see Sect. 4) un- 
til they either lose enough angular momentum to accrete onto the 
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star intact, or until they collide with a different fragment. Whether 
their orbits can fully contract before a catastrophic collision occurs, 
depends mainly on the size of the asteroid progenitor and the time 
required for circularisation. Larger asteroids produce more frag- 
ments and lead to faster collisions. Speed is key here and PR drag, 
the most suggested process for angular momentum loss, is clearly 
too slow. It was already shown by Veras et al. (2015b) that PR drag 
takes too long to accrete fragments above ~ 10 cm before the star 
cools down, let alone before they collide with other fragments. In 
our analysis of the tidal and material forces involved, we estimate 
that the upper limits to fragment sizes lie much higher, around 100 
m - 10 km depending mainly on the strength of the asteroid (see 
Sect. 2). 

There are, however, other processes that can contract orbits at 
a much greater pace than PR drag. One of these is the drag induced 
by fragment interactions with regions around the star that contain 
high concentrations of either gas or dust grains, for instance in the 
form of a compact disc that formed from in prior accretion event. In 
a recent work, Malamud et al. (2021) showed that drag at the frag- 
ment’s pericentre can even circularise km-sized bodies in several 
orbits, provided that a large second object already formed a mas- 
sive pre-existing disc around the star. If this disc only exists for a 
limited duration of time and is thick enough to circularise all but the 
very largest fragments, remaining km-sized boulders could survive 
into a new environment where they are relatively safe from colli- 
sional grind-down. Although its statistical significance has yet to be 
evaluated, this could form one channel to generate transiting mate- 
rial, as is observed around some systems (Vanderburg et al. 2015; 
Manser et al. 2019; Vanderbosch et al. 2020, 2021). The main ar- 
gument against the significance of this accretion channel is the lack 
of observations of gaseous emission lines or large IR excesses be- 
longing to the required pre-existing disc. 


7.5 Fork (5): A link between the dust circularisation speed 
and infrared excess 


In our suggested main channel, we continue with the collisional 
grind-down of fragments into dust. In this final phase, the speed of 
the dust circularisation leads to a split in possible observational out- 
comes. With our optically thin emission model, we show that slow 
dust circularisation in a sufficiently inclined disc via PR drag leads 
to detectable infrared excesses at higher accretion rates (> 107 g/s). 
This scenario could explain the minority of systems that show sig- 
nificant infrared excesses. If the inclination of the dust instead re- 
mains equal to the tiny value imparted during the tidal disruption 
event, the work done by stellar light becomes limited by the radial 
optical depth of the grains and the circularisation of the shaded dust 
grains slows down, ultimately limiting accretion onto the star and 
limiting the IR excess. 

However, the scenario of slow dust circularisation via PR drag 
cannot be used to explain the majority of systems that show no 
detectable infrared excess, even at high accretion rates. We sug- 
gest, therefore, that PR drag is not the only force that drives dust 
circularisation around most polluted white dwarfs. As was earlier 
suggested by Bonsor & Wyatt (2010), other drag forces - likely in- 
volving gas - are likely to play a key role. If the time required to 
circularise dust grains is reduced sufficiently by the gas drag, their 
accretion can occur without the accumulation of high grain abun- 
dances around the central star. In this manner, different circularisa- 
tion speeds of dust around different stars could break the propor- 
tionality between accretion rate and infrared excess. It is possible 
that the small quantities of gas required to accelerate the accretion 
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of dust grains are readily produced in the grind-down process itself. 
Indeed, Doppler tomography shows that some systems contain gas 
near the Roche radius, likely as a consequence of collisional pro- 
duction (Manser et al. 2016; Steele et al. 2021). 


7.6 Main road: model caveats and improvements 


Having discussed conceivable alternative paths to white dwarf pol- 
lution along with their physical selection criteria, we finally evalu- 
ate the model limitations of our suggested main road to accretion. 
The main uncertainty in the first stage relates to the fragment size 
distribution. As discussed in Sect. 2, the upper limit of the distribu- 
tion is limited by our uncertain knowledge of the material strength. 
At the lower end, it is not clear whether the smallest grains are in- 
deed produced in disruption events, and the slope is very poorly 
constrained. These factors severely limit the quantitative conclu- 
sions of our grind-down calculations, which are strongly related to 
our assumed size distribution. 

Despite the great uncertainties relating to the fragment sizes, 
we argue that it is worthwhile to evaluate the grind-down to ex- 
amine the process at order-of-magnitude scale and to investigate 
the trends and biases it involves. Our simple grind-down model of 
Sect. 5 should be interpreted in this way, rather than as an attempt 
to predict exact accretion rates. Firstly, it is based only on angu- 
lar differences induced by differential apsidal precession and does 
not include any other perturbing processes. Clearly these other per- 
turbing forces would play an important role. However, the preces- 
sion rates for differential apsidal precession are analytically known, 
such that these could be readily incorporated. We hypothesise that 
the general trends in collisional rates will follow a similar form. 
Secondly, the calculation only tracks catastrophic collisions and 
does not track the full collisional evolution of child orbits. This can 
be justified tentatively by the faster collisions of smaller fragments 
that result from the collisions but it remains an important limitation 
of the model. A more self-consistent evolution is possible to simu- 
late in theory but is numerically difficult considering the extreme 
eccentricity (~ 0.999) of the fragment orbits. Given our limited 
knowledge of the fragment size distribution, we did not consider 
that such a model would significantly improve our understanding 
of the processes involved. Although it could still be worthwhile to 
develop such a detailed model in the future, its predictive power 
will remain limited as long as the fragment size distribution after a 
tidal breakup event remains poorly constrained. 

Similarly, our calculation of the infrared excess in the dust ac- 
cretion stage is done in a simplified manner with the main goal to 
elucidate trends and show the two-dimensional morphology of the 
system rather than to predict precise excesses. Most importantly, 
we only performed calculations in the limiting case that the accret- 
ing dust is optically thin, whereas the discs become radially opti- 
cally thick if the rate of dust production is greater than ~ 107 g/s 
and the discs inclination remains small. Although it is clear that the 
inclination is directly linked to the observational outcome of dust 
accretion, further work is required to study how it evolves after the 
tidal disc is formed. 


8 SUMMARY AND CONCLUSIONS 


The main aim of this paper is to produce a road-map illustrating 
several potential routes for white dwarf pollution and to link these 
paths to observational outcomes (see Fig. 18). Our main route be- 
gins with the tidal disruption of a scattered asteroid and the for- 


mation of a highly eccentric tidal disc, followed by the collisional 
grind-down of fragments, which finally circularise and accrete onto 
the star as dust due to drag forces. Alternatives to this standard 
pathway for white dwarf pollution include a) the direct accretion of 
scattered asteroids/fragments of asteroids onto the white dwarf, b) 
the sublimation of dusty material and the accretion of gas, or c) the 
rapid circularisation of fragments via pre-existing compact discs. 
While accretion likely proceeds through a combination of these 
channels, the alternative a) is statistically very unlikely to occur, 
even with a planet re-scattering fragments of a disrupted asteroid. 
Channel b) will occur only for material scattered sufficiently close 
to the hottest white dwarfs, while c) occurs only following previous 
disruption events. 

Here we present detailed calculations of our suggested main 
road to white dwarf pollution. Our work includes simulations of 
collisional grind-down due to differential precession as well as a 
model for the infrared excess of the dust that is produced. Our main 
findings are that: 


(i) The size distribution of fragments in a tidal disruption event 
around a white dwarf can range as many as 10 orders of magnitude. 
The smallest bound fragments are no smaller than the limit set by 
radiation pressure at 0.1-10 um (Fig. 5), while fragments as large 
as 100 m-10 km also survive the tidal disruption depending on their 
material strength (Fig. 2). In the absence of a pre-existing compact 
disc or intense radiation, these larger fragments must be ground 
down before they can circularise and accrete by drag forces. 

(ii) Large asteroids produce more fragments when they disrupt, 
causing rapid collisional grind-down and generating short and in- 
tense bursts of dust production, whereas smaller asteroids grind 
down over longer periods of time. If subsequent dust accretion is 
fast, this biases observations to detect on-going accretion at inter- 
mediate rates by smaller asteroids (Fig. 13). Rare peaks in accretion 
rates from large asteroids are short-lasting and only probable to be 
detected in the atmospheres of DBZ stars with longer sinking time- 
scales (Fig. 12). 

(iii) Optically thin dust discs produce large amounts of infrared 
emission when their accretion rate exceeds 107 g/s. However, in or- 
der to remain completely optically thin at these high accretion rates, 
the inclination of the dust grains must be substantially increased 
beyond the tiny value imparted at the moment of tidal disruption. 
Infrared excesses at high accretion rates can be avoided by more 
rapid dust circularisation, for instance via enhanced drag due to the 
presence of gas near the disc's pericentre. 


DATA AVAILABILITY 


The simulation data that support the findings of this study are avail- 
able upon request from the corresponding author, Marc G. Brouw- 
ers. 
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